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Abstract 

We present a compositional compressible two-phase, liquid and gas, flow model 
for numerical simulations of hydrogen migration in deep geological radioactive 
waste repository. This model includes capillary effects and the gas diffusivity. The 
choice of the main variables in this model, Total or Dissolved Hydrogen Mass 
Concentration and Liquid Pressure, leads to a unique and consistent formulation 
of the gas phase appearance and disappearance. After introducing this model, we 
show computational evidences of its adequacy to simulate gas phase appearance 
and disappearance in different situations typical of underground radioactive waste 
repository. 

Keywords: Two-phase flow, compositional flow, porous medium, underground nuclear 
waste management 
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1 Introduction 

The simultaneous flow of immiscible fluids in porous media occurs in a wide vari- 
ety of applications. The most concentrated research in the field of multiphase flows 
over the past four decades has focused on unsaturated groundwater flows, and flows 
in underground petroleum reservoirs. Most recently, multiphase flows have gener- 
ated serious interest among engineers concerned with deep geological repository for 
radioactive waste. There is growing awareness that the effect of hydrogen gas gener- 
ation, due to anaerobic corrosion of the steel engineered barriers of radioactive waste 
packages(carbon steel overpacks and stainless steel envelopes), can affect all the func- 
tions allocated to the canisters or to the buffers and the backfill. The host rock safety 
function may even be threaten by overpressurisation leading to opening fractures in the 
host rock and inducing groundwater flow and transport of radionuclides outside of the 
waste site boundaries. 

Equations governing this type of flow in porous media are inherently nonlinear, and 
the geometries and material properties characterizing many situations in many appli- 
cations ( petroleum reservoir, gas storage, waste repository), can be quite irregular and 
contrasted. As a result of all these difficulties, numerical simulation often offers the 
only viable approach to modelling multiphase porous-media flows . In nuclear waste 
management, the migration of gas through the near field environment and the host rock, 
involves two components, water and pure hydrogen H2; and two phases "liquid" and 
"gas". Our ability to understand and predict underground gas migration is crucial to 
the design and to assessing the performance of reliable nuclear waste storages. This 
is a fairly new frontier in multiphase porous-media flows, and again the inherent com- 
plexity of the physics leads to governing equations for which the only practical way to 
produce solutions may be numerical simulation. 

This paper addresses one of the outstanding physical and mathematical problems 
in multiphase flow simulation: the appearance and disappearance of one of the phases, 
leading to the degeneracy of the equations satisfied by the saturation. In order to over- 
come this difficulty, we discuss a formulation based on variables which doesn't de- 
generate and hence could be used as an unique formulation for both situations, liquid 
saturated and unsaturated. We will demonstrate through four numerical tests, the abil- 
ity of this new formulation to actually cope with the appearance or/and disappearance 
of one phase in simple, typical but challenging situations, like the ones we met in un- 
derground radioactive waste repository simulations. 
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2 Modeling Physical Assumptions 



We consider herein a porous medium saturated with a fluid composed of 2 phases, 
liquid and gas, and two components. According to the application we have in mind, 
we consider the fluid as a mixture of two components: water (only liquid) and hydrogen 
(H2, mostly gas) or any gas with similar thermodynamical properties. In the following, 
for sake of simplicity we will call hydrogen the non-water component and use indices 
w and h for the water and the hydrogen components. 

We neglect the water vaporization since, in underground formations with high wa- 
ter pressure, the water vapor does not contribute significantly to the gas phase pressure 
. The water component is incompressible while the gas phase follow the ideal gas 
law. The whole fluid system is in thermal equilibrium and the porous medium is rigid, 
meaning that the porosity <t> is only a function of the space variable <t> = <t>(x); more- 
over, since hydrogen is highly diffusive we include the dissolved hydrogen diffusion in 
the liquid phase . 

The two phases are denoted by indices, I for liquid, and g for gas. Associated to 
each phase a € {l,g}, we have, in the porous medium, the phase pressures p a , the 
phase saturations S a , the phase mass densities p a and the phase volumetric flow rates 
q a . The phase volumetric flow rates are given by the Darcy-Muskat law: 

q, = -K(x)X,(S,) (V W - p /g ) , q g = -K(x)X g (S g ) (V Pg - p g g) , (1) 

where K(x) is the absolute permeability tensor, X a (S a ) is the a— phase relative mo- 
bility function, and g is the gravity acceleration; S a is the reduced a— phase saturation 
and then satisfies: 

S, + S g = l. (2) 
Pressures are connected through a given capillary pressure law: 

pc(Sg)=p g -pi- (3) 

From definition Q we notice that p c is a strictly increasing function of gas saturation, 
p' c (Sg) > 0, leading to a capillary constraint: 

Pg>Pl+Pc(0), (4) 

where p c {0) > is the capillary curve entry pressure ( see Figure]^. 

Since the liquid phase could be composed of water and dissolved hydrogen, we 
need to introduce the water mass concentration pj in the liquid phase, and the hydrogen 
mass concentration pj 1 in the liquid phase. Note that the upper index is the component 
index, and the lower one denotes the phase. We have, then 

Pi = Pi + Pi- (5) 

As said before, in the gas phase, we neglect the water vaporization and we use the ideal 
gas law: 

p g = C v p g , (6) 
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with C v = M h /(RT), where T is the temperature, R the universal gas constant and 
M h the hydrogen molar mass. Mass conservation for each component leads to the 
following differential equations: 

(Sipf) + div (pfqj = (7) 
(s lP >l + S gPg ) + div (p* q/ + p,qg + j* ) - ^* (8) 

where the phase flow velocities, q/ and q ? , are given by the Darcy-Muskat law (JTJ, & c 
are the source terms, and jf, c G are the c— component diffusive flux in liquid 



phase, as defined later in ( 13 i 



Assuming water incompressibility and that the liquid volume is independent of the 
dissolved hydrogen concentration, we may assume the water component concentration 
in the liquid phase to be constant, i.e.: 

P?=pf, (9) 

where p^. d is the standard water density. 

The assumption of hydrogen thermodynamical equilibrium in both phases leads 
to equal chemical potentials in each phase: Pg(T, Pg ,X^) — p^(T,p[,X^). Assuming 
that in the gas phase there is only the hydrogen component and no water, leads to 
Xg = 1; and then, from the above chemical potentials equality, we have a relationship 
p g = F(T,pi,X[ r ). Assuming that the liquid pressure influence could be neglected 
in the pressure range considered herein and using the hydrogen low solubility, p? <C 
P" = P'\l c '' we mav tr, en linearize the relationship between p g andX^, and finally obtain 
the Henry's law p g = K h Xf, where K is specific to the mixture water/hydrogen and 
depends only on the temperature T. Furthermore, using (|9]l and the hydrogen low 

, p''M w 

solubility, the molar fraction, X[\ reduces to p i,d M h ( se e eqs.(9)-(ll) in [2J) and the 
Henry law can be written as 

pi' = C h p gl (10) 

where C/, = HM h = p% d M h /(M w K h ); H is called the Henry law constant and is also 
depending only on the temperature. 

Remark 1 On the one hand the gas pressure obey the Capillary pressure law Q with 
the constraint Q, but on the other hand it should also satisfy the local thermodynam- 
ical equilibrium and obey the Henry law \10\ - More precisely if there are two phases, 
i.e. if the concentration, p^, is sufficiently high to have a gas phase appearance(S g > 0) 
, we have from and ([3]) : 

Pl =C h { Pl +p c (S g )). (11) 
Moreover, S g > with the constraint Q and the Henry's law ( 10 1, gives the constraint: 



tf>C H (pi+p e (0)). (12) 

But if the concentration, p^, is smaller than a certain concentration threshold ( see 
Figure^, then there is only the liquid phase (no gas phase, S g = 0), and none of all 



4 



S„ = 



pt<C h (p,+ Pr (0)) ^ 



> # / Is, > 

yy L? 1 

P9 = PI + Pc(Sg) 
>Pl+Po(0) 



Figure 1 : Phase diagram: Henry's law; localization of the liquid saturated S g — and 
unsaturated S g > states 



the relationships ([3]) or \12\ , connected to capillary equilibrium, applies anymore; we 
have only S g = 0, with pj< Qpg. 

There is then a concentration threshold line, corresponding to pj 1 = C/,(pi + p c {0)) in 
the phase diagram (Fig^, separating the one phase (liquid saturated) region from the 
two phase (liquid unsaturated) region. 

The existence of a concentration threshold line can also be written as an unilateral 
condition: 

0<S g <l, o<pj'<c hPg , s g (c hPg -p'n=o; 

which could be then used (see |4]) for designing a numerical scheme based on approx- 
imating a variational equation. □ 

The diffusive fluxes in the liquid phase are given by the Fick law applied to Z, w and 
to Xj 1 , the water component and the hydrogen component molar fractions (see eqs.(12) 
and (13) in [2 |). Using the same kind of approximation as in the Henry law, based on 
the hydrogen low solubility, we obtain, for the diffusive fluxes in this binary mixture 
(see Remark 2 and Remark 3 in |2|): 

jf = -<S>S l DVpj', jf = -j?, (13) 

where D is the hydrogen molecular diffusion coefficient in the liquid phase, corrected 
by the tortuosity of the porous medium. 

If both liquid and gas phases exist, (S g ^ 0), the porous media is said liquid unsat- 
urated and the transport model for the liquid-gas system can be now written as: 

*pF^ +div (pf 'q, - j*) = (14) 

®j t (S lP '! + C vPg S g ) + div (p/'q, + C vPg q g = (15) 

q / = -KA / (5 / )(v 7 , / -(pf + p/ 1 )g) , q g = -KX g (S g )(Vp g -C vPg g), (16) 

jj = -OS/DVp/'. (17) 
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But in the liquid saturated regions, where the gas phase doesn't appear, Si = 1 or 
S g = 0, the system ( 14 1-( 17 1 degenerates to: 

M 



div(pifq/-j* 



4>- 



dt 



-div(p, A q,+jf 



q, = -KA,(1) (V P; - (pf+pfig) , J? = -4>DVpf 



(18) 
(19) 



3 Liquid Saturated/Unsaturated state; a general for- 
mulation 

A typical choice for the two primary unknowns, in modeling immiscible two-phase 
flow, is the saturation and one of the phases pressure, for example S g and p/. But as 
seen above, in ( 14 1-( 19 1, this set of unknowns obviously cannot describe the flow in a 
liquid saturated region, where there is only one phase, and cannot take in account the 
gas dissolution since then the dissolved gas concentration, pj\ becomes an independent 
unknown. 



3.1 Modeling based on Total hydrogen concentration, p' t l ot 

To solve this problem, instead of using the gas saturation S g we have proposed, in 0, 
to use p/j„, the total hydrogen mass concentration, defined as: 



pl=S,p^+S g p h g . 



Defining 



with 



i(S g ) = C h (l-S g )+C v S g € [C h ,C v ]; 



a'(S g ) = C v — Ck = Ca > 0, 



(20) 
(21) 
(22) 



since C v > C%, from the assumption of weak solubility; we may then rewrite the total 
hydrogen mass concentration, p/j,,, defined in (20 1, as: 



h'tot 



a(S g )(pi+p c (S g )) ifS g >0 



Pi 



ifS g = 0. 



(23) 



With this new set of unknowns, p, ut and pi, the two systems of equations ( 14 1 



17) 



and ( 18 H 19 1 now reduce to a single system of equations: 



<S>pf ^ - div (pfKXi (S,) (V PI - (pf + p/')g) - OS/DVpf) = ^\ (24) 
- div (p* KA, (S,) (V PI - {pf + p*)g) 



- C vPg KX g (S g ) (Vp, + V Pc (S g ) - C vPg g) + *S,DVp} 



(25) 
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.QpM-JL-jjl-fa (A^V^+A^Vp^+BjKgJ -*P% d H dt ^ 



Now, if we want to study the mathematical properties of the operators in this system of 
equations, we should develop the above system of equations using S g = S g (pi,P( 0t ), 5/ = 
1 -S g = S,( Ph ft*,), and p* = pt{p h pi), t ), with 

djg = a ^^ 2]1 {P^>C;,(p/+p c (Q))} ds s = a ( S g) 1 {pA,>Q(p,+p f (Q))> (26) 

where 1^* >Ci < pi+p (q)u is the characteristic function of the set {p f * ( > C&(p/ + p c (0))}. 
As noted in section 2.5 in ||2), we have dS g /dpi < and dSg/dp^, > 0, when the gas 
phase is present. Then the system ( 14 i-([T5j) can be written : 

(27) 

-div(A 2,1 V^ + A 2 - 2 Vp*, +B 2 Kg) = J^. (28) 
Where the coefficients are defined by: 

A U (w,pL) =A/(5 g )p^K-4>(l-S g )Z)C /l iVI, (29) 
& l ' 2 (pupl) =-®(l-S g )^DC h I, (30) 

& 2A (PhPtlt) Hh(S g )pj l +X g (S g )C vPg N)K + <i>(l-S g )DC h NI, (31) 
A 2 - 2 (p/,P,t) =^(S,)^jq, ft K + *(l -S^^yDQI, (32) 

»i(ft,P&) =-A / (5,)pf [p;f + pf], (33) 
»2(Pi,/4) = - ik{S g )pt\pf +P 1 ,'} + Us s )clp 2 g ); (34) 

with I denoting the identity matrix and with the auxiliary functions 

PliPhPtot) = ™MChP g (phpL),Pwt), PgiPhPwt) = Pi+Pc(S g (j>i,pL))- (36) 



We should notice first that equation ( |28[ ) is uniformly parabolic in the presence 
of capillarity and diffusion; but if capillarity and diffusion are neglected, this same 
equation becomes a pure hyperbolic transport equation (see sec. 2.6 in [2 1). Then, if we 
sum equations (27 1 and (28 1 we obtain a uniformly parabolic/elliptic equation, which 
is parabolic in the unsaturated (two-phases) region and elliptic in the liquid saturated 
(one-phase) region. 

Remark 2 Simulations presented in sec. 3.2 in /|2y show that this last model with 
these variables, p^ ot , the total hydrogen mass concentration, and pi, the liquid phase 
pressure, could easily handle phase transitions (appearance and disappearance of the 
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gas phase) in two-phase partially miscible flows. However, in equations ( |27[ )-( |ZS| ), 
we should notice that both the coefficients A' J in operators, and the time derivative 
coefficients, can be discontinuous. For instance, only if the capillary pressure satisfies 



p' c (S g = 0) = +°o, as in the van Genuchten model, then all the coefficients in {27\-{28 \ 
are continuous; but if this condition is not satisfied they will be discontinuous. 

An other variant for the replacement of the saturation by p/j,,, is presented in 
[I], where relation (pO]iis written 

p& = (1 - S g )C hPg + S g C vPg , (37) 

and is then extended to both the two-phase and the one-phase region by making p g = p\ 
in the liquid saturated region (without the gas phase). This is leading, in the one-phase 
region where S g = and the Henry law does not apply, to extend the gas saturation by 



negative values (still defined by equation ( 37 1 as a function of the pressure and the total 



hydrogen concentration). After a necessary and ad hoc extension of the permeability 
and capillary pressure curves, out of the usual positive values of saturation, it is then 
possible to modeling both the one-phase flow and the two-phase flow with the same 
system of equations written with this extended saturation as main unknown, while using 
actually the total hydrogen concentration p/j,,. 

3.2 Modeling based on the hydrogen concentration in the liquid 
phase, pi 

We have seen that the variables pi and p/j„, introduced in the last section, can de- 
scribe simply the flow system, both in the one-phase and in the two-phase regions, 
independently of the presence of diffusion or capillary forces. But if we assume that 
the capillary forces are present we can choose an other change of variables in order 
to have a system of equations with continuous coefficients. Namely, using the inverse 
of the capillary pressure function, we may define the phase saturation as function of 
the hydrogen mass concentration in the liquid, p,' 7 , and of the liquid pressure, pf, and 
hence use them as main unknowns. With these two variables, pj 1 and pi the two sys- 
tems ([T4|-([T7]i and (p"8]l-([T9} are transformed in a single system of equations able to 
describe both liquid saturated and unsaturated flow. 

Since the capillary pressure curve S g H> p c (S g ) is a strictly increasing function we 
can define an inverse function / : R — > [0, 1], (see Fig. [2]), by 

m = {Pe\*) if*>fc(0) (3g) 

J y ' \ otherwise. v ' 

By definition of the function /, using (|T0]> and ( fT2) , we have: 

--S g , (39) 




and it is then possible to compute the gas saturations, S g , from pi and pj'. These two 
variables being well defined in both the one and two-phase regimes, we will now use 
them as principal unknowns. 



8 



Pc(Sg) 



i 4.... 



o 



^ Sg 











7T 



_ Pi 



Figure 2: Capillary pressure curve, p c = p g — p\, and inverse function 
Equations (fT4|i-(fT7|i with unknowns pi and p' ! can be written as: 



^5 ' 



-div(A 1 ' 1 Vp 



u Vpf+Bi 



where the coefficients are given by the following formulas: 

A 1 ' 1 = Xi(S g )p« d K, A 1 ' 2 = -Q(l-S g )DI 
A 2 ' 1 = A,6S,)P/*K, A 2 ' 2 = A A ,(S s ,)^p/'K + 4>(l 



B 1 = -WS g )p$ d (p$ d +p?) 



B: : -Ap^p/'.p-' +p,'<) A,(5 g )^(p/¥ 



v ^A\2 



with 



a* (S g ) 



a (S g 
Ch 



(-ft 



If we consider first, equation ( |4"Tj ), we may write it as 



ida*(S g ) dp, 



(40) 
(41) 

(42) 
(43) 

(44) 
(45) 

(46) 



' dp; 5? 

Moreover, from ( |4o*| i and because / and /' are positive, we have 



a*(S g )+p> 



da* (S g 
~dp~f 



and if the diffusion is not neglected, we have definite positiveness of the quadratic 
form A 2 ' 2 , in equation ( |4"Tj ); i.e. for any ^ 7^ 0, 

(A 2 ' 2 ^) =k g {S g )% P ^^+^{l-S g )D\^>0, 

and therefore equation ( |4"T| is strictly parabolic in p^. 
If we develop, equation (HOl as follows: 



®P% d f [ TT-Pl ] ^-div(A 1 ' 1 V P/ +A 1 - 2 Vp/ ! +B 1 Kg) 



C H " dt 




we have, for any ^ , 



and then positiveness of (A 1 ' 1 ^ ■§) and of (A 2,1 ^ ■ §). 
Moreover, 

^/'(^-w) >o. 



However, equations in system (|40]»-((4TJ are not uniformly parabolic/elliptic for the 
pressure pi, because the coefficients, A^^ A 21 , in front of Vpi in (|40"|i- ( |4"Tj ) tend to 



zero as S g — > 1 . 



Remark 3 It is worth noticing that this system ( |40| i-( |41| i, with variables p\ and p*, has 
interesting properties for numerical simulations in strongly heterogeneous porous me- 
dia. These two variables are continuous through interfaces separating different porous 
media with different rock types ( different absolute permeability, different capillary and 



permeability curves), as we will see in 4.3' which is absolutely not the case for the 
variables p\ and p^ ot . An other advantage is the continuity of all the coefficients A''f 
in ( |40| >-( |4T| > and the continuity of f in ( |41| l , even if p' c (S g = 0) = +°°. 



4 Numerical experiments 



In this last section, we present four numerical tests specially designed for illustrating 



the ability of the model described by equations ( 40 1-( 41 1 to deal with gas phase appear- 
ance and disappearance. Although all the computations were done using the variables, 
pi and p, , we are also displaying, for each test, the Saturation and Pressure level curves. 
These two last quantities are obtained after a post processing step using the Capillary 
Pressure law ([3]), equations (39 1, Henry's law (10 1, and following the constraints Q 
and ( 12 1 (see Figure[T]l. 
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The first test focuses on the gas phase appearance produced by injecting pure hy- 
drogen in a 2-D homogeneous porous domain Q. (see Figure[3]l, which is initially liquid 
saturated by pure water(water saturated). 

Because the main goal of all these numerical experiments is to test the model ef- 
ficiency, for describing the phase appearance or disappearance, the porous domain ge- 
ometry does not really matter and we will use a porous domain with a simple geometry. 
Consequently, we choose a simple, quasi- ID, porous domain (see Figure [4]) for the all 
three next tests . 

The test case number 2 is more complex, it shows local disappearance of the 
gas phase created by injecting pure hydrogen in a homogeneous unsaturated porous 
medium (initially both phases, liquid and gas, are present everywhere). 

The two last tests aim is to focus on the main challenges in simulating the flow 
crossing the engineered barriers, located around the waste packages. In the test case 
number 3, the porous medium domain is split in two parts with different and highly 
contrasted rock types, and like in the first one, the gas phase appearance is produced 
by injecting pure hydrogen in an initially water saturated porous domain. The test case 
number 4 addresses the evolution of the phases, from an initial phase disequilibrium to 
a stabilized stationary state, in a closed porous domain (no flux boundary conditions). 



Parameter 


Value 


e 


303 K 


d'i 


3 10~ 9 m 2 /s 




1 10~ 3 Pa.s 


Ih 


9 10- 6 Pa.s 


H(6=303K) 


7.65 10~ 6 mol/Pa/m 3 


M w 


10- 2 kg/mol 


M h 


2 10~ 3 kg/mol 


pf 


10 3 kg/m 3 



Table 1: Fluid parameters: phases and components characteristics. 

In all these four test cases, for simplicity, the porous medium is assumed to be isotropic, 
such that K = kl with k a positive scalar; and the source terms are assumed to be 
null: & w = and J^/, = 0. As usual in geohydrology, the van Genuchten-Mualem 
model for the capillary pressure law and the relative permeability functions are used in 





Mesh size range 


Time step range 


Test number 1 


2 m- 6 m W 


10 2 years - 5 10 4 years 


Test number 2 


1 m ( ** } 


10 2 years - 5 10 3 years 


Test number 3 


1 m ( ** ; 


10 2 years - 2 10 4 years 


Test number 4 


2 10~ 3 m (**' 


0.33 s- 16.7 10 3 s 



(*) Unstructured triangular mesh 
(**) Regular quadrangular mesh 



Table 2: Mesh sizes and time steps used in the different Numerical Test 
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imp 



r 



imp 



L, 



L, 



imp 



imp 



Figure 3: Test case number 1: Geometry a the 2-D porous domain, Q.. 



underground nuclear waste modeling, i.e. : 

' Pc = Pr (S;^ - l) ^ , h= - Vsie (l - (1 - S)l m T) 2 

and kg = — ^\-S l A\-S]l m ) 2m (47) 

Si — Si res 1 
with Sie = '— and m = 1 . 

Note that in the van Genuchten-Mualem model, we have no entry pressure, p c (0) = 0, 
but the presence of an entry pressure will not lead to any difficulty, neither from the 
mathematical point of view, nor for the numerical simulations. Concerning the other 
fluid characteristics, the values of the physical parameters specific to the phases (liquid 
and gas) and to the components (water and hydrogen) are given in Table [T] All the 
simulations, presented herein, were performed using the modular code Cast3m, Q. 
The differential equations system was first linearized by a quasi-Newton method and 
then discretized by a finite volume, implicit in time, scheme; with the discretization 
parameters (mesh size and time step) given in Table [2] 

4.1 Numerical Test number 1 

The geometry of this test case is given in Figure [3j and the related data are given in 
Table [3] A constant flux of hydrogen is imposed on the input boundary, T,,,, while 
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imp 



n 2 



Tout 



imp 



L, 



Figure 4: Test cases number 2, 3 and 4: Geometry of the quasi- ID porous domain, 



Dirichlet conditions pi — pi fiut , pf = are given on r oHf in order to have only the water 
component on this part of the boundary. The initial conditions, p\ = pi ou t and pf = 0, 
are uniform on all the domain, and correspond to a porous domain initially saturated 
with pure water. 



The main steps of the corresponding simulation are presented in Figure [5] 



Boundary conditions 
Initial condition 


Porous medium 


Others 


Param. 


Value 


Param. 


Value 


(j) K ■ v = on T imp 


k 


5 10- M in 1 


u 


200 m 


§ h ■ v = on r imp 




0.15 (-) 




20 m 


f ■ v = on r,„ 


Pr 


2 10 6 Pa 


Pl.out 


10 6 Pa 


(j) h ■ v = £ h on r,„ 


n 


1.49 (-) 


3 h 


9.28 mg/m^/year 


Pi = Pl.oui on r ot „ 




0.4 (-) 






p\ = on r out 


Sg,res 


o (-) 






Pl(t = 0)= Pl.out in^ 










pj t (t = 0) = 0mCl 











Table 3: Numerical Test case number 1: Boundary and Initial Conditions; porous 
medium characteristics and domain geometry; W and § h are denoting respectively 
the water and hydrogen flux. 

We observe in the beginning (see time t = 1200 years in Figure [5]) that all the 
injected hydrogen through T,-,, is totally dissolved in the liquid phase, the gas saturation 
stay null on all the domain (there is no gas phase). During that same period of time: 
the liquid pressure stay constant, the liquid phase does not flow, and the hydrogen is 
transported only by diffusion of the dissolved hydrogen in the liquid phase. 

Later on, the dissolved hydrogen accumulates around r,„ until the dissolved hy- 
drogen concentration pj 1 reaches the threshold pf — Cf,pi ( according to Figure [T^nd 
p c (0) = in[T]i, at time t = 1600 years, when the gas phase appears in the vicinity of 
r,„. Then this unsaturated region ( the two-phases, gas and liquid are present together) 
progressively expands and the liquid pressure, due to the compression by the gas phase, 
increases in the whole porous domain, causing the liquid phase to flow from r,„ to r„„ f . 
Consequently, after this time, t = 1600 years: the hydrogen is transported by convec- 
tion in the gas phase and the dissolved hydrogen is transported by both convection and 
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diffusion in the liquid phase. The liquid phase pressure increases globally in the whole 
domain until time t — 260 000 years (see Figure BJ, and it starts to decrease in the 
whole domain until reaching a uniform and stationary state at t = 10 6 years, in which 
the water component flux is null everywhere. 




pf at t = 1200 years p t at t = 1200 years S g at t — 1200 years 




p/ 1 at t = 4 10 4 years p t at t = 4 10 4 years S g at t = 4 10 4 years 




p* at / = 2 10 s years /?/ at f = 2 10 5 years 5 g at f = 2 10 5 years 




Pi at f = 10 6 years pi at t = 10 6 years at f = 10 6 years 



Figure 5: Numerical Test case number 1: Evolution of p/ ! ,the hydrogen concentration 
in the liquid phase; p\ the liquid phase pressure; and S g the gas saturation ; at times 
t = 1200,4 10 4 ,2 10 s and 10 6 years (from the top to the bottom). 
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4.2 Numerical Test number 2 



The geometry and the data of this numerical test are given in Figure |4] and Table [4] 
The porous medium is homogeneous and the initial conditions uniform; there is no 
need for defining two parts of the porous domain, £2; and i^; the parameter L\ will be 
considered as null. 

In this second test a constant flux of hydrogen is imposed on the input boundary 
Tin, while Dirichlet conditions p/ = pi ou t, p g — Pg.out are chosen, on Y out , such that 
Pl >ChPi, in order to keep the gas phase ( according to the phase diagram in Figure [IJ 
present on this part of the boundary. The initial conditions pi — pi l0U t and pf = C n pg tOUt 
are uniform and imply the presence of the gas phase ( S g > 0) in the whole domain. 



The main steps of the corresponding simulation are presented in Figures|6]and|7jivhere 
are presented the liquid pressure pi, the dissolved hydrogen molar density ( equal to 
pNM h ) and the gas saturation S g profiles at different times. 



Boundary conditions 
Initial condition 


Porous medium 


Others 


Param. 


Value 


Param. 


Value 


M '-v = Oon r imp 


k 


5 1(T 2U m 1 


L, 


200 m 


(j) h ■ v = on T imp 




0.15 (-) 


Ly 


20 m 


0"' • v = on r,„ 


Pr 


2 10 6 Pa 


Li 


m 


$ k - V = £ h on T in 


n 


1.49 (-) 


Pl,out 


10 6 Pa 


Pl = Pl.out on T oM 


Si, res 


0.4 (-) 


Pg.out 


1.1 10 6 Pa 


Pl — ChPg.out on Tom 


Sg.res 


o (-) 


£ h 


55.7 mg/m 2 /year 


Pi{t = 0) = pi. oul in £2 










pj>(t = 0)=C h p gtOUt in O. 











Table 4: Numerical Test case number 2: Boundary and Initial Conditions; porous 
medium characteristics and domain geometry. H and <j) h are denoting respectively 
the water and hydrogen flux. 



At the beginning, up to t < 1400 years, the two phases are present in the whole 
domain (see time t = 500 years on Figure [6]). The permanent injection of hydrogen 
increases both the two phase pressures and the gas saturation in the vicinity of r,„. The 
local gas saturation drop is due to the difference in mobilities between the two phases: 
the lower liquid mobility leads to a bigger liquid pressure increase, compared to the 
gas pressure increase; which is finally producing a capillary pressure drop (according 
to definition see Figure E}, and creating a water saturated zone. At time t = 1400 
years, the gas phase starts to disappear in some region of the porous domain (see time 
t = 1500 years, in Figure [7]l . 

Then, a saturated liquid region (S g = 0)will exist until time t = 17 000 years (see 
Figure [6]); and during this period of time, the saturated region is pushed by the injected 
Hydrogen, from Fi„ to T out . 

After the time t = 17 000 years, due to the Dirichlet conditions imposed on T out , 
the liquid saturated region disappears and all together the phases pressure and the gas 
saturation are growing in the whole domain (see the time t = 20 000 years in Figure|7]i. 

Finally the liquid pressure reaches its maximum at time t = 20 000 years and then 
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decreases in the whole domain (see the Figure[7]i. This is caused, like in the numerical 
test case number 1, by the evolution of the system towards a stationary state which is 
characterized by a zero water component flow. 





— 1= 


=0.00E+00 


— 1= 


=5.00E+02 


— 1= 


=1.50E+03 


t= 


=3.00E+03 


— 1= 


=8.00E+03 


— 1= 


=1.40E+04 


— t= 


=1.70E+04 



Figure 6: Test case number 2; L x = L 2 = 200 m: Time evolution of p h t (top right), the 
dissolved hydrogen molar density (pf/M ) (top left) and S g (bottom) profiles ; during 
the first time steps. 



4.3 Numerical Test number 3 

The geometry and the data of this numerical test are given in Figure |4] and Table [5] . 
Like in the Numerical Test number 2, a constant flux of hydrogen is imposed on the 
input boundary, r,„, while Dirichlet conditions pi — pi out , pf — are given on T out , in 
order to have only the liquid phase on this part of the boundary. The initial conditions, 
Pi — Pi, out an d Pi = 0, are uniform on all the domain, and correspond to a porous do- 
main initially saturated with pure water. Contrary to the two first numerical tests, the 
porous domain is non homogeneous, there are two different porous subdomains ill and 
H 2 ; L x = 200 m, L\ = 20 m and L 2 = 180 m. 

The simulation time of this test case is T = 10 6 years ;the discretization space mesh 
is 1 m; the time step is 10 2 years at the beginning and grows up to 2 • 10 4 years in the 
end of the simulation (see Tabl^2j. 
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50 100 150 200 50 100 150 200 

x (m) x (m) 



— t=1.70E+04 
— t=2.00E+04 
— t=3.00E+04 
— t=5.00E+04 
— t=1.00E+05 
t=2.00E+05 



"0 50 100 150 200 

X (m) 

Figure 7: Test case number 2; L x = Li = 200 m: Time evolution of the dissolved 
hydrogen molar density (p, A /M A ) (top right), p/(top left) and S g (bottom) profiles ; 
during the six last time steps. 




_ > s„ 



Figure 8: Saturation discontinuity at the interface of two materials with different cap- 
illary pressure curves; test case number 3. 
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Boundary conditions 


Porous medium 


Other 


initial condition 


Param. 


Value on £l\ 


Param. 


Value 


f ' • v = on T imp 




KT 18 


Lx 


200 m 


$ h • v = on Y imp 




0.3 (-) 


Ly 


20 m 


f" ■ V = on r,„ 




2 10 6 Pa 


u 


20 m 


(j) h ■ v = 3, h on Y in 


n 


1.54 (-) 


Pi, out 


10 6 Pa 


Pi = Pl.om on r OM , 


Sl,res 


0.01 (-) 


£ h 


5.57 mg/m^/year 


p/' = on T oW 


Sg.res 


o (-) 






Pl(t = 0)=pi, ou t on a 


Param. 


Value on O.2 






pf(f = 0) = 0on£2 


k 


5 10^ u m l 






4> 


0.15 (-) 








Pr 


15 10 6 Pa 








n 


1.49 (-) 








res 


0.4 (-) 








Sg.res 


(-) 







Table 5: Numerical Test case number 3: Boundary and Initial Conditions; porous 
medium characteristics and domain geometry. lv and ,! are denoting respectively 
the water and hydrogen flux. 



Figures[9]and[T0|represent the liquid pressure pi, the dissolved hydrogen molar density 
( equal to p, A /M") and the gas saturation S g profiles at different times. 

The main difference from the previous simulations (which were in a homogeneous 
porous domain) is the gas saturation discontinuity, staying on the porous domain inter- 
face x = 20 m; and due to the height of this saturation jump , we had to use a logarithm 
scale for presenting the gas saturation S g profiles . 

There are four main steps : 

• From to 3.8 • 10 4 years both the gas saturation and the liquid pressure stay 
constant in the whole domain while the hydrogen injection on the left side r,„ of 
the domain increases the hydrogen density level . 

• From 3.8 10 4 to5.4-10 4 years both the liquid pressure and the hydrogen density 
are increasing in the whole domain. The gas start to expanding from the left side 
of the domain r,„. The saturation front is moved towards the porous media 
discontinuity, at x — 20 m, which is reached at t =5.4- 10 4 years; see Figures[9] 
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the saturation front has 



From 5.4- 10 4 years to 1.3 • 10 5 years, see Figures 
crossed the medium discontinuity at x = 20 m and, from now, all the saturation 
profiles will have a discontinuity at x = 20 m. 



From 1.3 • 10 5 years to 10 6 years, see Figures 10 both the hydrogen density and 
the gas saturation keep growing while the liquid pressure decreases towards zero 
on the entire domain. The gas saturation front keeps moving to the right, pushed 



by the injected gas, up to x m 150 m at 10 years. 

Until the saturation front reaches the interface between the two porous media, for 
(f = 5.4 • 10 4 years), appearance and evolution of both the gas phase and the unsaturated 
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zone are identical to what was happening in the test case 1 (with a homogeneous porous 
domain) during the period of gas injection: the dissolved hydrogen is accumulating at 
the entrance until the liquid phase becomes saturated,at time (f > 3.8 10 4 years), letting 
the gas phase to appear. 

When the saturation front crosses the interface between the two porous subdomains 
(at x = 20 m and t — 5.4 ■ 10 4 years), the gas saturation is strictly positive on both sides 
of this interface and the capillary pressure curves being different on each side( see 
Table [5]) forces the saturation to be discontinuous for preserving the capillary pressure 
continuity on the interface. The capillary pressure continuity at the interface imposes to 
Pc > the Capillary Pressure in £li, and to pc, the Capillary Pressure in Q.2, to be equal 
on this interface, pc = P^P is satisfied only if there are two different saturations,on 
each interface side Sg l \ and S^: p^ {S^) = p^(S^) ; see Figure]^ 

In the same way as in the numerical test case number 1, the system tends to a 
stationary state 




I \ 

CO 



x (m) 




— 1= 


=0.00E+00 


— 1= 


=5.00E+03 


— 1= 


=2.00E+04 


t= 


=4.00E+04 


— 1= 


=4.60E+04 


— 1= 


=5.40E+04 



x(m) 



Figure 9: Test case number 3; L x = 200 m, L\ = 20m: Time evolution of the dissolved 
hydrogen molar density (p? /M h ) (top right), pi (top left) and S g (bottom) profiles ; 
during the first time steps. All the S g curves go to zero( although this cannot be seen 
using a logarithmic scale) . 
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— 1= 


=5.40E+04 


— t= 


=7.00E+04 


— t= 


=1.00E+05 


t= 


=1.30E+05 


— t= 


=2.00E+05 


— 1= 


=5.00E+05 


— t= 


=1.00E+06 



Figure 10: Test case number 3; L x = 200 m, L\ = 20m: Time evolution of the dissolved 
hydrogen molar density (p' ! /M /, )(top right), p\ (top left) and S g (bottom) profiles ; 
during the last seven time steps. All the S g curves go to zero( although this cannot be 
seen using a logarithmic scale) . 
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4.4 Numerical Test number 4 



This last numerical test is different from all the precedent ones; it intends to be a simpli- 
fied representation of what happens when an unsaturated porous block is placed within 
a water saturated porous structure. The challenge is then: how the mechanical balance 
will be restored in a homogeneous porous domain , which was initially out of equilib- 
rium, i.e. with a jump in the initial phase pressures? 



The initial liquid pressure is the same in the entire porous domain; £2 ,pi\ = pi 2 , 
and in the subdomain £2i the initial condition, (p/.i = p g .\ in Table|6j corresponds to a 
liquid fully saturated state with a hydrogen concentration reaching the gas appearance 
concentration threshold (p g = p\ and p^ 1 = Chp g , see Figure [IV In the subdomain Q.2 
the initial condition(p/ 2 ^ p ? ,2and p g 2 ^ p g ,i Corresponds to a non saturated state 
(see Table[6|. 

The porous block initial state is said out of equilibrium, because: 
if this initial state was in equilibrium , in the two subdomains £2i and Q.2, the local 
mechanical balance would have made the pressures, of both the liquid and the gas 
phase, continuous in the entire domain £2 . 

For simplicity, we assume the porous medium domain £2 is homogeneous and all 
the porous medium characteristics are the same in the two subdomains £2i and £22, and 
corresponding to concrete. 

The system is then expected to evolve from this initial out of equilibrium state 
towards a stationary state. 

We should notice that, in order to see appearing the final stationary state, in a rea- 
sonable period of time, we have shortened the domain £2 ( L x = lm), taken the porous 
media characteristics , and set the final time of this simulation 7/,„ at 7/;„ — 10 6 s « 
1 1.6 ) days. The complete set of data of this test case is given in Table [6] 



Boundary conditions 
initial condition 


Porous medium 


Other 


Param. 


Value 


Param. 


Value 


W • v = on d£l 


k 


1(T 1S m z 


Lx 


1 m 


<j) h - V = on da 




0.3 (-) 


Ly 


0.1 m 


Pi{t = Q)=pi.\ onQ] 


Pr 


2 10 6 Pa 


U 


0.5 m 


pf(t = 0)=C hPgA on Cli 


n 


1.54 (-) 


Pl,l 


10 6 Pa 


Pl(t = 0)= pi, 2 on Q.2 


$l,res 


0.01 (-) 


PgA 


10 6 Pa 


pf(f = 0) =C hPg .2 on Q.2 


Sg.res 


(-) 


Pl,2 


10 6 Pa 








Pg,2 


2.5 10 6 Pa 



Table 6: Data of the numerical test number 4 : boundary and initial conditions ;domain 
geometry. The porous medium domain £2 is homogeneous, all the porous medium 
parameters are the same in the two subdomains £2; and £22; <j) w and 0' 1 are denoting 
respectively the water and hydrogen flux . 

The space discretization step was taken constant equal to 2 ■ 10~ 3 m and the time 
step was variable, going from 0.33 s in the beginning of the simulation to 16.7 • 10 3 s at 
the end of the simulation (see Table |2j. Figures 1 1 and[T2"]represent the liquid pressure 
pi, the dissolved hydrogen molar density ( equal to p^/M' 7 ) and the gas saturation S g 
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0.2 0.4 0.6 0.8 1 

x (m) 



Figure 12: Numerical test case number 4, L x = 1 m, L\ = 0.5 m: Time evolution of the 
the dissolved hydrogen molar density (Pf /M h ) (top right), p/ (top left) and S g (bottom) 
profiles; during the five last time steps. 
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profiles at different times. 

There are essentially two steps: 



For < t < 1 .92 • 10 5 s (see Figure Hi, the initial gas saturation jump moves from 
x = 0.5 m, at t — and reaches T,,,, the left domain boundary, at t = 1.92 • 10 5 s. 
During this movement, the saturation jump height (initially w 0.16 ) decreases, 
until approximately 0.03, when it reaches the left boundary r,„. In front of this 
discontinuity there is a liquid saturated zone, S g = 0, and in this zone both the 
liquid pressure and the hydrogen density are spatially uniform (see Figure [TT[ 
top ). But, while the hydrogen density remains constant and equal to its initial 
value, the liquid pressure becomes immediately continuous and starts growing 
quickly (for instance, p/(t = 10 3 s) sa 1.6 • 10 6 Pa), and then more slowly until 
t = 1.3 ■ 10 5 s, when it starts to slightly decrease. 



In Figure 1 1 located on the gas saturation discontinuity, there are both a high 
contrast in the dissolved hydrogen concentration (this concentration stays how- 
ever continuous, but with a strong gradient, as seen in the top right of Figure 
and a discontinuity in the liquid pressure gradient (see the top left of Figure 



For 1 .92 • 10 s s < t < 10 6 s = 7y,„ (see Figure 12 1, all the entire domain is now 
unsaturated (S s >0). The liquid pressure, the hydrogen density and the gas satu- 
ration profiles are all strictly monotonous and continuous, going towards a spa- 



tially uniform distribution, corresponding to the stationary state(see Figure 12 1. 



As expected, the system initially out of equilibrium (discontinuity of the gas pres- 
sure), becomes immediately again in equilibrium (the gas pressure is continuous)and 
evolves towards a uniform stationary state (due to the no mass inflow and outflow 
boundary conditions). Although the liquid pressure and the dissolved hydrogen den- 
sity are immediately again continuous for t > , the hydrogen density still have a 
locally very strong gradient until t — 1 .92 • 10 s s. 



At first, and at the very begining(sa 10 2 s), see top left of Figurelll] only the liquid 
pressure evolves in the liquid saturated zone. Due to a gas pressure in the unsaturated 
zone higher than in the liquid saturated zone (S g = 0; p H — 2.5MPa > p/ = IMPa, for 
the initial state in Tabl^6]l, and due to the no flow condition imposed on r,„, the liquid in 
the saturated zone is compressed by the gas from the unsaturated zone. Then, a liquid 
gradient pressure appears around the saturation front and makes the liquid to flow from 
the liquid saturated zone towards the unsaturated one, and then the gas saturation front 
to move in the opposite direction. 

The very strong hydrogen density gradient (until t = 1.92 • 10 5 s), located on the 
saturation front, is due to the competition between the diffusion and the convective flux 
of the dissolved hydrogen around the saturation front: the water flow convecting the 
dissolved hydrogen, from left to right, cancels the smoothing effect of the gas diffu- 
sion propagation in the opposite direction. On the one hand the diffusion is supposed 
to reduce the hydrogen concentration contrast, by creating a flux going from strong 
concentrations (in the unsaturated zone) towards the low concentrations (in the liquid 
saturated zone), and on the other hand the flow of the liquid phase goes in the opposite 
direction (left to right, from S g = to S g > 0). Once the disequilibrium has disap- 
peared, the system tends to reach a uniform stationary state determined by the mass 
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conservation of each component present in the initial state (the system is isolated, with 
no flow on any of the boundaries). 



5 Concluding remarks 

From balance equations, constitutive relations and equations of state, assuming ther- 
modynamical equilibrium, we have derived a model for describing underground gas 
migration in water saturated or unsaturated porous media, including diffusion of com- 
ponents in phases and capillary effects. In the second part, we have presented a group 
of numerical test cases synthesizing the main challenges concerning gas migration in 
a deep geological repository. These numerical simulations, are based on simplified but 
typical situations in underground nuclear waste management; they show evidence of 
the model ability to describe the gas (hydrogen) migration, and to treat the difficult 
problem of correctly following the saturated and unsaturated regions created by the gas 
generation. 
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